Molecular Dynamics Simulation of Shear Moduli for Coulomb Crystals 
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Torsional (shear) oscillations of neutron stars may have been observed in quasiperiodic oscilla- 
tions of Magnetar Giant Flares. The frequencies of these modes depend on the shear modulus of 
neutron star crust. We calculate the shear modulus of Coulomb crystals from molecular dynamics 
simulations. We find that electron screening reduces the shear modulus by about 10% compared 
to previous Ogata et al. results. Our MD simulations can be extended to calculate the effects of 
impurities and or polycrystalline structures on the shear modulus. 
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Recently quasi-periodic oscillations (QPOs) have been 
observed in the tails of Magentar Giant Flares [I] [2]. 
These flares are extremely energetic 7-ray bursts from 
very strongly magnetized neutron stars. The QPOs have 
been interpreted as shear oscillations of the crust 
If this interpretation is correct, the QPO frequencies 
could provide detailed information on neutron stars and 
their crusts [5j . The frequencies of shear modes depend 
on shear moduli of neutron star crust which is a Coulomb 
solid. Ogata et al. [S] have calculated shear moduli us- 
ing Monte Carlo simulations. In this paper wc improve 
on the Ogata et al. results by presenting molecular dy- 
namics simulations with much better statistics and we 
include the effects of electron screening. 

Our results for shear moduli are preparation for later 
molecular dynamics calculations of the breaking strain of 
neutron star crust [7]. This breaking strain determines 
the maximum height of neutron star "mountains" be- 
fore they collapse under their own weight. Mountains 
on rapidly rotating neutron stars may efficiently radiate 
gravitational waves |1] . These waves could limit the spin 
frequencies of accreting stars and may be detectable with 
large scale interferometers 0. In addition the breaking 
strain may be important for crust breaking models of 
Magnetar Giant Flares (lOj . 

In the crust of a neutron star electrons form a very 
degenerate relativistic gas. The ions are completely pres- 
sure ionized and have Coulomb interactions that are 
screened at large distances by the slightly polarizable 
electron gas. The interaction potential between two ions, 
v{r), is assumed to be |1F, 



with rie the electron density. The total potential energy 
is Vfot = X]i<j "^(^ij)- Charge neutrality ensures that 
Ue — Zn where n is the ion density. The ions are assumed 
to form a classical one component plasma (OCP) that can 
be characterized by the Coulomb parameter F, 



F = 



(3) 



This parameter is the ratio of a typical Coulomb to ther- 
mal energy and the ion sphere radius a — [3/(47rn)]^/^ 
characterizes the separation between ions. The OCP is 
expected to freeze for F > 175. 

To calculate shear moduli, we follow the procedure of 
Ogata et al. [S]. The change in free energy with defor- 
mation 5F can be expressed in terms of elastic constants 
cii, C12 and C44, 



= ^(Cll - Cl2)u% 



+ CiiUikUki (i^ k) , (4) 



and Uife describes the strain. 

Under a deformation, the coordinates of an ion get 
mapped to r^. 



(5) 



fc=i 



We consider six deformations Di {i — 1...6) that conserve 
the volume to order e^. 



Uxx = e + 



(6) 



v{r) 



(1) 



where the ions have chage Z, r is the distance between 
them, and the electron screening length Ae is 



rl/2 



2e(37r2ne)i/3 
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(2) 



D2 ■■ Uyy ^ e + -e , Uxx = Uzz = 
3 2 

D^'- = e + -e , Uxx = Uyy = 
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Dq : U^_^ = = - , Uyy^— (11) 

For each deformation Dm we calculate a corresponding 
expectation value /„i (m — 1...6), 

where V is the system volume. At zero temperature, this 
reduces to = {(fVtotlde^)/V. 

For a body centered cubic crystal one has [6], 

h^h^h^ 36n = 3(cn - C12) (13) 

and 

/4 = /5 = /e = C44. (14) 

Here cn, C12, and C44 are elastic constants. In practice 
we calculate all six independently and average to de- 
termine 611 and C44. The angle averaged shear modulus 
is IS], 

Mcff = (26ii+3c44)/5. (15) 

If neutron star crust involves many crystal domains of 
random orientation, then /iog is the appropriate elastic 
constant to determine the speed of shear waves. 

The shear modulus is sensitive to the very long range 
tails of the interactions. To study this we cut off the 
potential at a large distance i?cut i 

v{r) ^ i;cut(r) = [v{r) - i)(i?cut)]e(i?cut - r). (16) 

We have subtracted a constant so that Vcutir) is con- 
tinuous at r = -Rcut- In Fig- [l] we plot the elastic 
constants &11 and C44 versus i?cut for a perfect bcc lat- 
tice at zero temperature. This figure was calculated as- 
suming Z = 29.4. We note that the ratio of Ae to a 
is Ae/a = 5.41/Z^/'^ independent of density. We see 
that one must go to very large i?cut > 12Ae to calcu- 
late both 611 and C44 accurately. For i?cut ^ 00 we have 
^off = 0.1108(nZ2e2/a). This is 8% smaller than the 
value /Zoff = 0.1194(nZ^e^/a) that Ogata et al. [5], cal- 
culate in the limit Ag ^ 00. We conclude that electron 
screening , neglected in ref. f6j, reduces /iog by about 
10%. 

We now describe our MD simulations at finite tempera- 
tures. For simplicity we work at a density n = 7.18 x 10"^ 
fm^'' and Z = 29.4. Our results can be scaled to other 
densities at a given value of F. Our results can also be 
approximately scaled to other values of Z, at fixed F. 
This is because, although the ratio Ag/a changes with Z, 
this change in screening has only a small effect on the 
shear modulus. We evolve the system with the velocity 
Verlet algorithm [T3] using a time step 5t — 25 fm/c. 
Starting from T = and a perfect bcc lattice we increase 
the temperature to T = 0.1 MeV and evolve the sys- 
tem for typically 100000 MD steps (2.5 x 10^ fm/c) to 
reach thermal equilibrium. Next we evolve for a further 




FIG. 1: (Color on line) Elastic constants 3feii and C44 versus 
cutoff distance i?cut for a perfect bcc lattice at zero tempera- 
ture. The cutoff distance is in units of the electron screening 
length Ae. 



250000 MD steps (6.25 x 10^ fm/c) storing configurations 
for later calculations of elastic constants. The tempera- 
ture is then raised by of order 0.1 MeV and the process 
repeated. We keep the system at a fixed temperature 
(approximately) by periodically rescaling the velocities. 
These MD simulations are done in an undistorted cubic 
box using periodic boundary conditions. 

We calculate fm by averaging over 1000 configurations, 
each separated by 250 MD steps (6250 fm/c). To mini- 
mize finite size effects we calculate T4ot by summing over 
all 27 nearest periodic images. Thus ion i is assumed to 
interact not only with ion j at its original position but 
also with 26 more images of j where the x, y, and z coor- 
dinates are independently shifted by 0, +1, or — /, with I 
the box size. The derivatives in Eq. [12] are approximated 
using a five point numerical formula. We note that the 
MD trajectories have been calculated using periodic dis- 
tances (involving only the single nearest periodic image 
of a given ion) to save time, while the derivatives have 
been calculated by summing over 27 images to minimize 
finite size effects. 

Table |T] presents results for simulations using N = 3456 
ions and no cutoff i?cut = 00. Statistical errors only are 
indicated in parentheses. We caution that 611 may have 
significant errors from finite size and other systematic ef- 
fects. Indeed Fig. [T] suggests that finite size effects could 
be large for this small system. However 611 only makes a 
small contribution to /Xoff- Therefore /ies in Table]!] may 
be more accurate. We fit the values of /Ltes in Table^with 
a simple analytic formula that is valid for all F > 175, 

98 7 72p2 
Meff « (0.1106 -—)(n ). (17) 
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TABLE I: Shear Moduli for MD simulations with N = 3456 
ions. 

r bii (nZ^e^/a) C44 (nZ^e^/a) figs {nZ'^e^/a) 



00 0.0220 
834 0.0209(2) 
417 0.0194(2) 
278 0.0202(4) 
200 0.0154(5) 
175 0.0158(8) 



0.1699 

0.1617(3) 

0.1517(3) 

0.1410(5) 

0.1253(10) 

0.1152(10) 



0.1107 

0.1054(2) 

0.0988(2) 

0.0927(3) 

0.0813(6) 

0.0755(6) 



TABLE IL Shear Moduli for MD simulations with iV = 9826 
ions using a cutoff -Rcut = 13.9Ae. 



00 0.0212 
834 0.0208(2) 
200 0.0177(5) 



0.1700 

0.1602(2) 

0.1224(6) 



0.1105 

0.1045(1) 

0.0805(4) 



This fit has an error < 2%. 

To study finite size effects we have performed addi- 
tional simulations with larger systems. Table [IT] presents 
results for simulations with N = 9826 ions using a cut- 
off i?cut = 13.9Ae. For this larger system and for finite 
r, /icff is about 1% smaller in Table |n] than in Table |T] 
Therefore we estimate finite size effects in Table |TT] to be 
of order 1%. Figure [2] plots these results for /Xeff and 
in addition shows results for very small simulations with 
N — 1024 ions, where finite size effects are large. Finally, 
Fig. |2]also shows the Monte Carlo results of Ogata et al. 
[3]. These results are about 10% larger than our results 
at large F and have much larger statistical errors. 

Ogata et al. neglect electron screening ^ 00. 
At zero temperature we have performed calculations for 
larger values of Ag and extrapolated to Ae ^ 00. Note 
that we can not directly calculate for Ag = cxd. Our ex- 
trapolated results are consistent with Ogata et al. There- 
fore we conclude that electron screening reduces ficS by 
about 10%. The speed of shear waves is proportional to 
the square root of ficS- Therefore electron screening re- 
duces the shear speed by about 5%. This will slightly 
lower the frequency of torsional oscillations of neutron 
star crusts. 

In future work, we will study the impact of impurities 
on /Xeff by explicitly including them in our MD simula- 
tions |12j . We expect impurities to lower the shear modu- 
lus because they reduce the uniformity of the crystal. We 
will also study the effect of polycrystalline structure on 
flag with larger scale MD simulations that include multi- 



ple crystal domains. These multiple domains could also 
lead to a lower effective shear modulus. Finally, we will 
calculate the breaking strain by slowly deforming the sim- 
ulation volume and calculating the resulting stress. The 
breaking strain is important for the maximum height of 
mountains on neutron stars that could be important for 
gravitational wave radiation. In addition, the breaking 
strain is important for star "quakes" that may trigger 
Magnetar Giant Flares. 

In conclusion, we have calculated the shear modu- 
lus of a Coulomb plasma using MD simulations. The 
shear modulus is important for the frequencies of tor- 
sional oscillations of neutron star crusts. Our results 
for the angle averaged shear modulus ficS are, iicS ~ 
(0.1106 - 28.7/ri-3)(nZ2eVa)- Here n is the ion den- 
sity, Z the ion charge, and a the ion sphere radius, 
a — (3/47m)^/^. This formula is accurate to about 2% 
and valid for Coulomb parameter F > 175. Our results 
are about 10% smaller than Ogata et al because we in- 
clude electron screening. 
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FIG. 2: (Color on line) Angle averaged shear modulus ficB 
versus Coulomb parameter F for MD simulations involving 
TV = 1024, 3456, and 9826 ions. Also shown are Monte Carlo 
results from Ogata et al. ^ that omit electron screening. 
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